Heat conduction in a three dimensional anharmonic crystal 
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We perform nonequilibrium simulations of heat conduction in a three dimensional anharmonic 
lattice. By studying slabs of length TV and width W , we examine the cross-over from one-dimensional 
to three dimensional behavior of the thermal conductivity We find that for large A^, the cross-over 
takes place at a small value of the aspect ratio W/N . ^From our numerical data we conclude that 
the three dimensional system has a finite non-diverging «: and thus provide the first verification of 
Fourier's law in a system without pinning. 

PACS numbers: 



Macroscopic behavior of heat transport in the hnear 
response regime is governed by Fourier's law 



J 



-kVT(x), 



(1) 



where J, VT are respectively the heat current density 
and temperature gradient at the position x, and k, is the 
thermal conductivity. This implies diffusive behavior of 
heat. What are the necessary and sufficient conditions 
for the validity of Fourier's law ? This question is a 
longstanding unsolved problem [1]. For solids one starts 
with the description in terms of a harmonic crystal where 
heat conduction takes place through lattice vibrations or 
phonons. Scattering of the phonons can occur due to 
phonon-phonon interactions (ie anharmonicity in the 
interactions) or by impurities (e.g isotopic disorder, de- 
fects) f^. For one dimensional systems, from a large 
number of numerical and analytical studies it is now es- 
tablished that these scattering mechanisms are insuffi- 
cient in ensuring normal diffusive transport. Instead one 
finds anomalous transport [3, 4], one of the main signa- 
tures of this being that the thermal conductivity k, in 
such systems is no longer an intrinsic material property 
but depends on the linear size N of the system. A power 
law dependence k ~ A^" is typically observed. For two 
dimensional anharmonic crystals a k ~ In(A^) divergence 
of the conductivity is predicted from various analytical 
theoriesjl, |6| and also from an exactly solved stochastic 
model [7|, but the numerical evidence for this so far is 
inconclusive [1, 0| • A recent experiment has reported the 
breakdown of Fourier's law in nanotubes [lo| while an- 
other experiment on graphene flakes ^ij also indicates a 
divergence of k. 

For systems with pinning {i.e an external substrate 
potential) and anharmonicity, Fourier's law has been ver- 
ified in simulations on one and two dimensional systems 
Q. There is a strong belief that Fourier's law should 
be valid in three dimensional_^(3i?) systems, even with- 
out pinning. A recent work examined heat transport 
in a iD disordered harmonic crystal. Analytical argu- 
ments showed that heat conduction in the system was 



sensitive to boundary conditions. For generic boundary 
conditions a finite conductivity was predicted but this 
could be numerically verified only for the pinned case. In 
this letter we investigate the effect of anharmonicity on 
heat conduction in ordered crystals. Through extensive 
simulations of a 3Z3 anharmonic crystal we give strong 
numerical evidence for normal transport and the validity 
of Fourier's law in this system. 

Model. — We consider a 31? cubic crystal with a scalar 
displacement field a;„ defined on each lattice site n = 
(ni,n2,n3) where ni = 1,2, ...,A^ and 112 — — 
1, 2, W. The Hamiltonian is taken to be of the Fermi- 
Pasta-Ulam (FPU) form: 



n n,e 
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where e denotes unit vectors in the three directions. We 
have set the values of all masses and harmonic spring 
constants to one and the anharmonicity parameter is 1/. 
Two of the faces of the crystal, namely those at ni = 1 
and ni = N, are coupled to white noise Langevin type 
heat baths so that the equations of motion of the particles 
are given by: 

in = -^[{Xn- Xn+e) + v{Xn - a^n+e)^ ] 
e 

+ 5m,i{-ix^ + ri^) + 5m,N{-lx^ + ri^) ■ (3) 

The noise terms at different sites are uncorrelated 
while at a given site the noise strength is specified by 
(77^'«(t)77^.«(t')) = 2jTL,RSit - t') , where n and Tr 
are the temperatures of the left and right baths and we 
have chosen units where the Boltzmann constant kg = 1. 
Fixed boundary conditions were used for the particles 
connected to the baths and periodic boundary condi- 
tions were imposed in all the other directions. We sim- 
ulate these equations using a velocity- Verlet algorithm 
[isll and calculate the heat current and the temperature 
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profile in the nonequilibrium steady state of the crystal. 
The heat current jn from the lattice site n to n + ei 
where ei = (1,0,0), is given by j„ = (/„,„+§! in+ei), 
with /n,n+ei being the force on the particle at site n + ei 
due to the particle at site n. In our simulations we cal- 
culate the average current per bond given by 



J 



N-l W 

ni — 1 712 1^3 — 1 



We also calculate the average temperature across layers 
in the slab and this is given by T„^ = {1/W'^) I]«2,«3 *n- 
Simulation details. — In all our simulations we set 
f = 2 and = 2, Tr = 1. We first address the question 
of the dependence of J on the width W of the system 
and the nature of the cross-over from ID behaviour, for 
small values of the ratio r = W/N, to true 3D behavior 
for W/N ~ 1. The numerical results are given in Fig. ([T}. 
We see that for any fixed length N, the value of J de- 
creases as we increase W but saturates quickly to the 3D 
value. The cross-over width Wc is seen to increase slowly 
with A^. The inset shows that as we increase A'^, the cross- 
over from ID to 3D behavior takes place at decreasing 
values of r and presumably in the thermodynamic limit 
N ^ oo, the cross-over occurs at r — > 0. Thus our study 
suggests that Wc ~ N"" with < a < 1. A similar result 
was obtained by Grassberger and Yang for a 2D FPU 
system. 

Next we look at the dependence of J on N for the 
3-D case. The fast cross-over from ID to 3D behaviour 
implies that we can extrapolate the results for small r 
to estimate the true value of the 3D current (at r = 1). 
Thus we can get results for quite large values of N from 
simulations on systems with small widths. For sizes up 
to iV = 128 we obtained data for W = N. For the largest 
system size, namely N = 16384 we have data for W ~ 
16. We show our results for the N dependence of k in 
Fig. ([2]). There are three sources of error in the values of 
current: (i) numerical errors, arising from the finite time 
discretization value {dt — 0.001), and from rounding off 
errors; (ii) statistical errors arising from averaging over a 
finite number of time steps; and (iii) errors arising from 
the extrapolation of the small aspect ratio (r) results to 
the 3D case. The error from (iii) was taken to be the 
difference in current values for the two largest widths 
studied. For smaller system sizes we verified that the 
numerical error was much smaller than the statistical and 
extrapolation errors and we assume that this is true also 
at larger system sizes. The error-bar for each data point 
plotted in Fig. 1^ is the larger of errors from (ii) and 
(iii). 

The slope of the n versus N curve is decreasing slowly 
with TV and a straight line fit to the last three points gives 
an exponent a ~ 0.09 ± 0.01. For comparison we also 
show in Fig. ^ the ID and 2D data for the FPU system. 
The 2D results are from data ioi N x N samples for 
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FIG. 1: Plot of the heat current J versus width W for different 
fixed values of the length A*'. The inset plots J versus the 
aspect ratio r — W/N. 



100 



10 



1 

- QJV 


.\d'. 

Mi:--.. -3D\ 

° iS s * 

-B * t + * 


1 1 1 1 1 1 1. 

I 

z 

z 


3).l 








z N £ 


. - ' 

• - ' ' • ' ID ; 


- 1 




. 2D ■ 







, , , , ■ , 3D ," 



28 

N 



2^2 2" 2 



16 



FIG. 2: Plot of K-versus-A' in different dimensions. The inset 
shows the running slope on = din k/ din N as a function of 
A''. The dashed line is a guide to the eyes. 



systems up to A^ = 2048 while for larger sizes the results 
shown are extrapolated values from small width data. In 
ID we get a w 0.33 [Ii] while in 2D we get a sa 0.22. 
In the inset of Fig. ([2]) we have plotted the running slope 
defined as ajv = din k/ din N against system size. From 
this we see that while the slopes in ID and 2D tend 
to saturate, the 3D slope seems to be decreasing. The 
3D slope can be fitted by the dashed line with a power 
law form. This suggests that the asymptotic system size 
behaviour will give a = implying diffusive transport 
and validity of Fourier's law. 

One of the remarkable features of ID systems with 
anomalous heat transport is the form of the steady state 
temperature profile obtained in these systems. Typically 
one finds that the temperature profile is concave upwards 
in part of the system and concave downwards elsewhere 
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and this is true even for small temperature differences 



14l . Il5j . This means that the temperature gradient is 



non-monotonic as a function of distance across the sam- 
ple. In Fig.® we plot the temperature profiles for the 
ID, 2D and 3D samples. We see that the variation of 
the temperature gradients are non-monotonic in both ID 
and 2D while in 3D they are monotonic. The inset in 
Fig. (l3])shows that the 3D temperature profile is concave 
upward everywhere. We have also confirmed that the 
profile becomes more linear on decreasing the tempera- 
ture difference between Tl and Tr. This again supports 
our finding based on the size-dependence of the current, 
that heat transport in 31? is diffusive while in lower di- 
mensions it is anomalous. 

Finally we look at the temperature dependence of 
thermal conductivity. Temperature and nonlinearity are 
highly correlated [16|, and temperature dependence can 
be understood from the nonlinearity dependence of ther- 
mal conductivity. We note that Eq. ([3]) leads to the 
scaling relation sJ{T, AT, siy) — J{sT,sAT,v), where 
T ^ {Tl + Tr)/2, at = Tl - Tr, and s is an ar- 
bitrary scale factor. Taking the limit AT -> 0, this 
gives the scaling relation for thermal conductivity as 
k{T,si>) — k{sT,v). Putting v = 1 and s ~ v,we then 
get 



k(T, v) = k{i>T, 1). 



(4) 



Thus the thermal conductivity is a function of vT. One 
may expect that large ly suppresses heat currents due 
to enhancement of phonon-phonon interactions. Hence 
from the scaling (jH) we expect that k must also decrease 
with increasing T. To check this, we show the dependence 
of the heat current on i/T for a 32 x 32 x 128 system with 
a small temperature difference AT = 0.1. In Fig.Q, we 
compared two cases: one with v = 2.0 fixed and T var- 
ied, and another with T — 1.0 fixed and h' varied. We 
find that current decreases as a function of vT, consistent 
with the scaling relation Eq. ([4]). We note that Fourier's 
law (P leads to (fT/dx^ = - J"^ k'^ {dn / dT) and so the 
decrease of k in the region T g [1.0,2.0] with v — 2.0 
is consistent with the concave curve in the 3D temper- 
ature profiles. Interestingly at large anharmonicity the 
current does not seem to go to zero but instead appears 
to saturate to a constant value. At low temperatures the 
effect of anharmonicity becomes weaker and we expect 
the conductivity to increase, eventually diverging in the 
limit T — )■ 0. It is difficult to numerically access the low 
temperature regime since the mean free path becomes 
large and one would need much larger system sizes to see 
diffusive behaviour. 

Summary and Discussion. — In summary, we have 
given the first numerical evidence for the validity of 
Fourier's law of heat conduction in an anharmonic crys- 
tal in three dimensions. This confirms the belief that 
in three dimensions anharmonicity is a sufficient condi- 
tion for normal transport. This is not a necessary con- 



dition since, for example, a 3D pinned disordered purely 
harmonic crystal also shows normal transport |12l |. Our 
conclusion was based on three evidences. The first is the 
system-size dependence of the thermal conductivity, the 
second is temperature profile, and the third is the consis- 
tency between temperature profile and temperature de- 
pendence of conductivity. It has been known that the 
one-dimensional FPU system shows slow convergence of 
the thermal conductivity to it's asymptotic behavior [l3 |. 
Here we show that this is also the case in 3D. Unlike ID 
and 2D, the running slope of the size dependence of n 
in 3D showed decreasing behavior even at the largest 
system size and this gives us a clear signature for finite 
K. The temperature profiles in 3D are completely differ- 
ent type from the ID and 2D case where nonmonotonic 
behavior of the gradient is robust even for small tem- 
perature differences. We note that a recent simulation 
of heat conduction in the 3D FPU crystals reported di- 
verging thermal conductivity (the reported exponent is 



about 0.221) [17|. The reasons for this is probably be- 



cause of the small values of anharmonicity used in those 
simulations and also the much smaller system sizes that 
were studied (maximum size in that study was N = 256). 
In 2D we find a divergence of the conductivity with an 
exponent a ~ 0.22 which is similar to the value obtained 
in i. 

For a sample of fixed length N we find that the cur- 
rent density decreases on increasing its width W and the 
cross-over from ID to 3D behaviour takes place at a value 
Wc ~ iV" with < a < 1. This has implications for ex- 
periments measuring thermal conductivity of nanowires 



20[. If the cross-over width were independent of N 
and the width of the nanowire larger than it, then the 
thermal conductivity of long nanowires could well be fi- 
nite and not diverge as expected for true ID systems. 
On the other hand, since the cross-over width gradually 
increases with increasing N, a gradual transition from 
3D-like to ID behavior will take place when the cross- 
over width is comparable to the width of nanowires. This 
scenario is an interesting system size effect that may be 
observed in experiments on nanowires. 
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FIG. 3: Plot of temperature profiles for a ID system, a N x N 
2D system and a x x A'^ 3D system with different aspect 
ratios r — W/N. Temperature profiles for the three aspect 
ratios overlap with each other. The inset shows that the 3D 
temperature profile is concave upward everywhere. 
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FIG. 4: Demonstration of scahng property (H)) for 32 x 32 x 128 
system with AT = 0.1. Thermal conductivities decrease as 
increasing temperature T or nonlinearity v. This tempera- 
ture dependence of k explains the slightly concave curve in 
temperature profiles in 3-D (see Fig.(|3}). 



